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Abstract: We present an implementation of the next-to-leading order hadronic produc- 
tion of a boson in association with a pair of massive bottom quarks in the framework 
of POWHEG, a method to consistently interface NLO QCD calculations with shower Monte 
Carlo generators. The process has been implemented using the POWHEG BOX, an automated 
computer code that systematically applies the POWHEG method to NLO QCD calculations. 
Spin correlations in the decay of the boson into leptons have been taken into ac- 
count using standard approximated techniques. We present phenomenological results for 
Wbb — )• Ivbb production, at both the Tevatron and the LHC, obtained by showering the 
POWHEG results with PYTHIA and HERWIG, and we discuss the outputs of the two different 
shower Monte Carlo programs. 
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1. Introduction 

The production of a boson in association with a pair of massive bottom quarks {b and 
b), contributing to both the W + lb-jet and W + 2 6-jcts signatures, represents both an 
interesting Standard Model signal and one of the most important background processes for 
single-top production and Higgs searches. 

The cross sections for W boson production with bottom quarks has been measured at 
the Tevatron pp collider at Fermilab by both the CDF [1] and DO [2] Collaborations. As 
more data will be collected and analyzed by the Tevatron Collaborations, we will gain in- 
creasing precision in the cross section measurements and we will have a unique opportunity 
to test and improve the theoretical description of heavy-quark jets at hadron colliders by 
performing a thorough comparison between the Tevatron experimental data and existing 
theoretical predictions. Studying the same cross sections in the very different kinematic 
regimes available at the LHC pp collider will then be of great interest and will represent a 
crucial test of our understanding of QCD at high-energy colliders. 

Moreover, the production of a 14^ boson with one or two b jets represents a crucial 
irreducible background for both single-top production (j>p,pp — )■ tb,tb) and Higgs-boson 
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associated production (WH), followed by the decay H — t- bb. We remind that WH asso- 
ciated production is the most sensitive Higgs-boson production channel at the Tevatron, 
while it is a difficult, but very important, channel at the LHC, where, in particular kine- 
matic regions (boosted H), can provide essential additional signal for the detection of a 
low-mass Higgs boson [3, 4, 5]. 

It is therefore crucial to have the W + 6-jets background theoretically under good 
control. Several steps have already been taken towards this goal. At the parton level, the 
production oi a W boson with up to two jets, one of which is a 6 jet, has been calculated 
including next-to-leading-order (NLO) QCD corrections in the variable-flavor scheme [6], 
while the production of a 1^ boson with two b jets has been computed at NLO in QCD 
using the fixed-flavor scheme, first in the massless 6-quark approximation [7, 8, 9, 10] and 
more recently including full 6-quark mass effects [11, 12, 13, 14]. The two calculations 
have been combined in ref. [15] to provide the most accurate theoretical predictions for 
11/^-1-1 6-jet production. The comparison with the experimental measurement of the total 
cross section for W plus at least one b jet shows a clear discrepancy [1, 16, 17] which 
should be further investigated, given the importance of this background, in particular for 
Higgs boson searches. A few predictions for W + 2b jets at NLO in QCD can be found in 
ref. [11] for the Tevatron and in ref. [13] for the LHC, and are available to the experimental 
community for comparison. 

The Standard Model prediction for Wbb production could be further improved by 
properly interfacing the parton level NLO calculation with a parton shower (PS) simula- 
tion. An event generator of this nature should also be beneficial in understanding other 
experimental systematics for which parton shower simulations are relied upon, including 
the transverse momentum and rapidity distributions of b jets and hard non-6 jets, as well 
as the angular separation and invariant-mass distribution of the two 6-jet system, often 
used in experimental analyses to enhance the signal to background ratio. 

In recent times, the construction of these NLO-I-PS event generators has become viable: 
see, e.g., MCONLO [18] and POWHEG [19, 20]. The effectiveness of the POWHEG approach has 
been demonstrated successfully and studied in some detail through its application to a 
substantial array of hadron collider processes (see, for example, [21, 22, 23, 24, 25, 26, 27, 
28, 29, 30, 31, 32, 33, 34, 35]). In this paper we interface the NLO calculation of the cross 
section for Wbb production to a shower Monte Carlo program within the POWHEG framework. 
This is the first time that such calculation has been performed: more specifically, we have 
implemented this process using the POWHEG BOX [36] , a general computer code framework for 
embedding NLO calculations into shower Monte Carlo programs according to the POWHEG 
method. Spin correlations in the decay of the vector boson into leptons have been taken into 
account using standard approximated techniques [37]. We have checked that we have very 
good agreement between our approximated result and the NLO calculation for Wbb lubb 
production of ref. [14], where spin correlations have been taken into account exactly. 

The paper is organized as follows. In section 2 we give a brief description of the POWHEG 
implementation of the Wbb process and of the needed ingredients. Since we use the POWHEG 
BOX to implement our process, we refer the reader to the POWHEG BOX publication [36], and 
we report here only those aspects of the implementation that are particularly relevant to 
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the process in question. In section 3, after briefly recalling how POWHEG generates an event 
and how the POWHEG BOX deals with the presence of Born-zero configurations, we illustrate 
how we have implemented the decay of the W boson. In section 4 we present and discuss a 
few results for the Tevatron and the LHC, where the POWHEG-generated events are showered 
by PYTHIA and HERWIG. Finally, in section 5, we summarize our results and give our final 
remarks. 

2. Construction of the POWHEG implementation 
2.1 The next-to- leading order cross sections 

The NLO QCD corrections to qq' — > Wbb production consist of both one-loop virtual 
corrections to the tree level processes depicted in fig. 1 and one-parton real radiation 
from both the initial- and final-state quarks, i.e. qq' — t- Wbb + g. At the same order, the 
QdiQQ') ~^ Wbb + q'{q) processes also need to be included. The O (a^) virtual corrections 




Figure 1: The tree- level Feynman diagrams for qq' — > Wbb. 

consist of self-energy, vertex, box and pentagon diagrams. Dimensional rcgularization is 
used to regularize both UV and IR singularities. The UV singularities are cancelled by 
introducing a suitable series of counterterms. We renormalize the wave functions of the 
external massive quarks in the on-shell scheme, and the strong coupling constant ag in 
the mixed scheme of ref. [38], in which the heavy-flavour top loop is subtracted at zero 
momentum, while all the other flavours arc subtracted in the MS scheme. Self-energy, 
vertex, box and pentagon diagrams contain IR divergences that combine and cancel against 
the soft and collinear divergences in the real-emission corrections, when computing infrared- 
safe quantities with renormalized parton distribution functions (pdf). More details on the 
calculation of the O (a^) virtual corrections can be found in refs. [11, 12, 13]. 

While in refs. [11, 12, 13] the W boson is not decayed and is produced on-shell, the 
POWHEG BOX implementation of Wbb production includes the leptonic decay of the W^, in an 
approximated way, as described in section 3.5. To implement this, the analytic results of 
refs. [11, 12, 13] have been modified to lift the W on-shell condition and the invariant mass 
of the W boson has been generated according to a Breit-Wigner distribution function. 
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2.2 The POWHEG BOX ingredients 

We have implemented the NLO QCD eross seetion for Wbb production into the POWHEG 
BOX environment by providing the following ingredients: 



i. 
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trie 


list of all flavour structures of the Born processes; 


ii. 


the 


list of all flavour structures of the real processes; 


iii. 


the 


squared Born amplitude B and the colour correlated ones Bij^; 


iv. 


the 


Born phase space; 


V. 


the 


squared real matrix elements for all relevant partonic processes; 


vi. 


the 


finite part of the virtual corrections computed in dimensional regularization: 


vii. 


the 


Born colour structures in the limit of a large number of colours. 



The symmetric colour correlated Bij amplitudes, according to the particle labelling used 
in the POWHEG BOX (g(l) q'{2) W{3) 6(4) 5(5)), are given by 



for the case where qisa quark. For the case where the first incoming particle is an antiquark, 
we have the same results with Bu and B15 exchanged. Here Ca = 3 and Cp = 4/3 are the 
Casimir invariants of the colour SU(3) representation. We have taken the squared Born 
amplitude B from ref. [11] while the squared real amplitudes have been generated using 

MadGraph [39]. 

The assignment of colour flow for the two Feynman diagrams at the Born level is 
straightforward and unambiguous, and follows directly the propagation of quarks. 

2.3 Validation of the NLO code 

In the POWHEG BOX framework, it is possible to compute NLO distributions, taking advan- 
tage of the fact that the POWHEG BOX computes automatically all the counterterms needed 
to regularize the real distributions. It is then possible to check that, in the collinear and 
soft limits, the real amplitude has the correct behaviour, and that it is consistent with the 
Born cross section and the colour- and spin-correlated Born amplitudes. 

The NLO distributions obtained within the POWHEG BOX, using a variant of the FKS 
subtraction scheme [40], have been checked against the codes developed in refs. [11, 12, 
13] which use instead a phase-space slicing method with two cutoffs, to extract soft and 
collinear singularities analytically [41, 42, 43], as well as independent calculations of both 
virtual and real corrections, performed with several tools based on FORM [44], TRACER [45] 
and MAPLE codes. Full agreement has been found in all the studied distributions, assuring 
us that the ingredients provided to the POWHEG BOX are correct and consistent. 

^We notice that the B'^" spin-correlated Born amphtudes are zero, since there are no external gluons at 
the Born level. 
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3. Theoretical introduction 



In the POWHEG formalism, the generation of the hardest emission is performed first, at full 
NLO accuracy, and subsequent radiation is generated using shower Monte Carlo programs 
(see ref. [20] for a more detailed description of the method). In the following, we will briefly 
summarize a few features of the POWHEG method that will be useful in view of the discussion 
of sec. 3.4. 

3.1 Generation of the underlying Born configuration 

The first step of the generation process is the construction of the underlying Born kine- 
matics, i.e. the generation of the Born momenta, distributed according to the function 

Bi^n) = 5(*n) + V{^n) + J C^^rad R{^n+l) , (3.1) 

where B, V and R are the Born, virtual and real contributions to the NLO cross section, 
specifies the kinematics of the underlying Born event with n final-state particles, $rad 
the radiation variables and the real-emission variables *n+i = {^m *^'rad} are parametrized 
in terms of the underlying-Born and radiation ones. 

3.2 Generation of the radiation variables 

Once the momenta of the underlying Born have been generated, the POWHEG method pro- 
ceeds to generate the radiation, i.e. $rad) starting from the POWHEG cross section for the 
generation of the hardest emission 

da = S(*„) A + A (*„, (*„+i)) :^^^^d$rad| , (3.2) 

where values of A;^ ($n+i) < Pt™ not allowed (here ~ 1 GeV, i.e. of the order of 
a typical hadronic scale) . The Sudakov form factor A is given by 



A(*„,PT)=exp|- J d^.ad ^(fex (^n+l) - Pt)| . 



(3.3) 



The function k-j. (^n+i) should be equal, near the singular limit, to the transverse momen- 
tum of the emitted parton relative to the emitting one. The cross section in eq. (3.2) has 
the following properties: 

• at large k-j. it coincides with the NLO cross section up to next-to-next-to-leading order 
terms. In fact, in the large transverse-momentum region, eq. (3.2) can be written as 

da = 5(*„) ^^^^ d^n rf^rad , (3.4) 

since the Sudakov form factor approaches 1 in this region. This differs from the pure 
NLO result because of the presence of the factor 

= (3.5, 
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For processes that get large radiative corrections, the 0{as) term can be in fact of 
order 1, giving then a harder spectrum for the generated radiation. 



• It reproduces correctly the value of infrared-safe observables at NLO. Thus, also its 
integral around the small A;t region has NLO accuracy. 

• At small k-T it behaves no worse than standard shower Monte Carlo generators. 

3.3 Tuning of the real contribution 

In POWHEG it is possible to tune the contribution to the real cross section that is treated 
with the Monte Carlo shower technique. This was pointed out first in ref. [19], where the 
POWHEG method was formulated, and then implemented in the POWHEG BOX as a general 
feature. In this way, the enhancement in eq. (3.5) can be controlled, if necessary. In fact, 
the real cross section can be split into two positive contributions 

R = Rs + Rf, (3.6) 

such that Rf has no singularities (soft or collinear) and only Rg is singular in the infrared 
regions. In previous implementations [22] the separation was achieved using a function F 
of the transverse momentum of the radiation, ^ F ^ 1, that approaches 1 when the 
transverse momentum of the radiated parton vanishes, and such that 

Rs = RF, (3.7) 
Rf = R[l-F]. (3.8) 

The generation of the radiation is then done by POWHEG using only the divergent contribu- 
tion Rs and eqs. (3.2) and (3.3) become 

da = Bsi^n) Ci*n| A, (*n,Pr") + A. (*n, (*n+l)) ^^^^^^^^jj 

+ i?/(*„+l)d*„d$rad, (3.9) 
Bsi^n) = 5(*n) + Vi^n) + J d^^ad ^.(*n+l) , (3.10) 

A,(*„,pT) = exp|-y" rf^rad ^ff^'^) ^ ^(^T (^n+l) -Pt)} . (3.11) 

The contribution Rf, being finite and positive, is generated with standard NLO techniques, 
and fed into a shower Monte Carlo as is. 

3.4 Wbb production 

As for Higgs-boson production in gluon fusion [22], also in Wbb production the NLO 
corrections are very large, independently of the choice of the renormalization and factor- 
ization scale [11, 13]. For example, using a renormalization and factorization scale equal 
to ^ = niw + 2m;, and input parameters as specified in sec. 4, for W^bb production at the 
LHC with ^/s = 14 TeV, the LO cross section corresponding to the distributions shown in 
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rub [GeV] 


LO [pb] 


NLO [pb] 


K 


0.1 


140.6 ±0.3 


412 ±8 


2.9 


1.0 


135.5 ± 0.2 


381 ±2 


2.8 


10 


37.11 ±0.03 


98.2 ± 0.3 


2.6 


100 


0.5240 ± 0.0003 


0.961 ± 0.003 


1.8 



Table 1: Values of the LO and NLO cross sections, as well as their respective K factors, for W~bb 
production at the LHC with 14 TeV, using = mw + 2m;, as renormalization and factorization 
scale and varying the quark mass nib by three orders of magnitude. 

fig. 2 is 81.32 pb, wliile the NLO cross section is 222.88 pb, with a K factor (ratio of NLO 
over LO total cross section) of 2.74, while at the Tevatron, the cross sections are 10.43 pb 
at LO, and 19.84 pb at NLO, with a K factor of 1.9. 

We can attribute these large NLO corrections to at least two reasons: the opening 
of a new gluon-initiated channel in the real contributions at NLO (e.g. qg — t- Wbb + q), 
that is likely to be more important at the LHC than at the Tevatron [11, 12, 13], and the 
presence of large logarithms of the mass of the bottom quark, related to final-state coUinear 
singularities for massless quarks, now regularized by the quark mass. To illustrate this, we 
have collected in table 1 the K factors for several values of the mass of the final-state heavy 
quark, for the LHC at 14 TeV, and using renormalization and factorization scales equal 
to fj, = mw + 2m;,. One can see that the K factors are large and, when the quark mass 
increases from 0.1 GeV to 100 GeV, the K factor decreases from 2.9 to 1.8. We expect the 
large collinear logarithms to affect distributions even more than total cross sections. 

Before presenting a few results for Wbb production, we would like to address a further 
issue related to this process: Wbb production has a Born zero when a relativistic W is 
emitted parallel with respect to the incoming beam, and consequently the gluon (which 
undergoes the splitting into a bb pair) is emitted in the opposite direction. Along the 
incoming beam, angular-momentum conservation cannot be preserved, and this configura- 
tion is kinematically suppressed. This was in fact the case for W production, where there 
is a zero in the Born cross section if the outgoing lepton from W decay is anti-parallel to 
the incoming quark [21]. In fact, due to the left-handed nature of the W boson coupling, 
we have a violation of angular-momentum conservation along the incoming beam. 

Despite the fact that -B(*°) can vanish in some kinematic regions this kinematics 
can be generated since the Bs{^n) term in eq. (3.10) can be different from zero, due to 
the real term. Prom eq. (3.9), we see that, in this case, away from the Sudakov region 
(i.e. where ~ 1), the real contribution may be enhanced by a factor Bg/B, now large 
since B — > 0. The POWHEG BOX has a built-in mechanism to deal with this problem: is 
chosen to vanish in the regions where R differs too much (by more than a factor of 5) from 
its collinear or soft approximation, which are proportional to the underlying Born cross 
section. The contribution Rf = R — Rg, being non-singular, is then added independently. 
To activate this mechanism, the POWHEG BOX flag bornzerodamp has to be set to 1. This 
must then necessarily be done for the process we are studying. 
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Figure 2: Transverse-momentum distributions for the W boson, , and the hardest radiated 
non-fe jet, p;^, in W^bb production at NLO in QCD, for both the LHC with ^/s — 14 TeV (upper 
plots) and the Tevatron (lower plots). The dotted blue lines represent the results of the pure NLO 
QCD calculation. The POWHEG results obtained with the mechanism to protect from Born zero 
not activated and with no separation of the hard-radiation region are shown in dashed black lines. 
In solid red, the results with the mechanism to protect from Born zero activated and with the 
separation of the hard radiation collinear to the bottom quarks, as described in sec. 3.4. 



In fig. 2, we have plotted the transverse momentum of the W boson and of the hardest 
radiated jet at the Tevatron and at the LHC. Jets are reconstructed using the anti-fex 
algorithm [46] with R = 0.4, and jets are recombined using the default recombination 
E scheme [47]. No other cuts are applied to the events. The values of all the coupling 
constants, masses and physical parameters for the generation of these events can be found in 
sec. 4. In dashed black lines, the POWHEG hardest-emission results when the bornzerodamp 
flag is set to false. It is evident the effect of the enhancing factor Bg/B with respect to 
the corresponding NLO results (dotted blue lines). 

In addition to the activation of the mechanism to protect from the Born zeros, we have 
decided to separate out the region of hard-gluon emission collinear to one of the final-state 
massive quarks from the part of the real contribution that is treated by the Monte Carlo 
shower techniques (i.e. generated through the Sudakov form factor), and to handle it with 
standard NLO techniques, as described in sec. 3.3. In fact, this region is responsible for the 



-8- 



presence of large logarithms of the mass of the quark and we do not want these logarithmic 
terms to get further enhanced by the B/B ratio of eq. (3.5). We would like to stress that, 
strictly speaking, this region is not singular, since the mass of the quark regularizes it, but 
it would be a singular region if the mass of the quark were exactly zero. 

In order to separate out the region of hard-gluon emission collinear to the final-state h 
quarks, we have chosen the following form for the function F of eqs. (3.7) and (3.8) 

(i/d)^ + (iM)^ + (i/d-,)- ^ 

where 

d = E"^ {I - cos^ 9) , (3.13) 

and where k is the momentum of the radiated gluon with energy E, forming an angle 6 
with the positive direction of the incoming beam, and an angle Of, with the outgoing b 
quark, in the centcr-of-mass frame. The momentum of the b quark is ki,, with energy Ef, 
and three-momentum k^. The d^ term is similar to the d^ one except for exchanging b with 
b in all the kinematic variables. We have set c = 1 in the code, but higher values can be 
used too. 

The function F in eq. (3.12) has the following properties: 

1. it approaches 1 in the singular region, i.e. when the emitted gluon is parallel to the 
incoming beams or soft (d — > 0) , assuring that the singular region is treated with the 
Monte Carlo shower technique; 

2. it becomes small when the radiated parton is hard and collinear to the b or the b 
quark. In fact, when the radiated gluon is hard and collinear to one of the two heavy 
quarks, the db or di terms reach their minimum value. 

The distributions obtained with the mechanism to protect from the Born zeros and with 
the F function of eq. (3.12) are plotted as solid red lines in fig. 2. We find very good 
agreement with the NLO curves, at least in the region where we expect this to happen. 
The expected disagreement in the low-px jet region is due to the fact that the NLO curve is 
divergent in this region, while the resummed POWHEG result feels the effect of the Sudakov 
form factor and goes correctly to zero. In addition, we have verified that the major role 
in getting this agreement is played by the activation of the mechanism to protect from 
Born-zero configurations, while the effect of the separation of the region of hard gluons 
parallel to the massive b quarks plays only a minor role. 

3.5 boson decay 

Since we have used the analytic calculation of the virtual corrections of refs. [11, 12, 13], 
that treats the W boson as stable, we are not in a position to have all the spin correlations in 
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the leptonic M^-boson decay products correctly accounted for. We can instead use standard 
techniques to implement the decay in an approximated way [37]. In this approximation, 
spin correlations are not accurate to NLO in the whole phase space, but are correct to 
NLO for hard real emissions and to leading order in the soft and collinear region. In order 
to achieve this, we have produced a W boson with invariant mass M distributed according 
to the Breit-Wigner function 

2^ mw^W (rt-\K\ 

where niw and Tw are the pole mass and width of the W boson. Using the POWHEG method, 
a Born-like, or real-like event, is generated with an undccaycd W boson, whose invariant 
mass is M, and whose kinematics is parametrized by a set of variables that we call ^>u) 
where "u" stands for "undecayed" . 

The procedure that we are going to follow can be easily illustrated if we recall that the 
squared matrix elements are connected to the concept of probability. We then rephrase 
the procedure that we have used in a probabilistic language. The differential probability 
distribution of the decay variables is proportional to (we neglect the overall normalization 
factor that ensures that the integral of the differential probability distribution is 1) 

dP($d) - Mdi^d) d^d = Md(*u, ^w^iu) d^w^iu , (3.16) 

where Aid is the squared amplitude corresponding to the decayed process, with finite-width 
effects fully taken into account^. For consistency, the squared amplitude Aid must include 
only resonant diagrams (i.e. diagrams where the W momentum equals the sum of the / and 
ly momenta). In writing eq. (3.16), we have parametrized the kinematics of the process 
for the production and decay of the W boson in terms of the undecayed variables $u and 
of a set of variables describing the W decay, ^w^w- Equation (3.16) implicitly defines 
d^W^l^. Similarly, the differential probability distribution of the undecayed variables is 
proportional to A^u(^u) '^*^'u- The problem is then to determine the probability distribution 
of the variables that parametrize the decay, ^w^lui given the probability distributions for 
$d and To solve this problem we use the fact that the joint probability of two events 
A and B can be written in terms of the conditional probability 

P{AuB) = P{A\B)xP{B), (3.17) 

that in our case becomes {A corresponds to the generation of the decay variables for W ^ lu 
and B corresponds to the generation of the undecayed Whh event) 

Addi^d) d^d ^ dPi^w^iu I ^'u) X 7Wu(^u) d<Pn , (3.18) 



^The decayed tree- level squared amplitudes Add have been obtained using MadGraph [39], both for the 
Born and for the radiative processes. 
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where dP(^\Y^ii, \ <I>u) is the infinitesimal probabihty distribution of the variables ^w^lu^ 
given the kinematics of an undecaycd process We can then write, using eq. (3.16), 

dP{^w-.i. I $u) ^ = AU^^^ 



d^w^iv ■ (3.19) 



dP{^w~^iu I *^'u)/'^^Ty->-ii/ is the distribution function we are looking for. To generate effi- 
ciently ^w^lv, distributed according to (3.19), we use the hit-and-miss technique and so 
we need to find an upper bound for the ratio M.a{^\i-,^w-^Iv) / We have used as 
upper bounding function the expression 

Ud{M\ ^w^i^) = iVd 2, '2x^2 Mw^UM^) , (3.20) 

where is a normalization factor, Mw^lu is the squared decay amplitude corresponding 
to the W ^ Iv decay and mi is the charged-lepton mass. One can predict the appropriate 
value for the normalization factor or compute it by sampling the decay phase space 
(^w-^lv and comparing with the exact expression, in such a way that the inequality 

holds. The veto algorithm is then applied as follows: 

1. first one generates a point in the phase space 

2. then a random number r in the range [0, Ud{M'^, (^w^lu)] is generated; 

3. finally, if r < A^d(^u5 ^w^lv) / M.\i{^\i)i the kinematics of the decay is kept and the 
event is generated. Otherwise the algorithm goes back to step 1. 

At the end of this procedure, the kinematics ^w^lv of the W decay is generated, and, 
together with the POWHEG- generated undecayed variables the kinematics of Whh event 
followed by the decay of the W boson becomes available. 

In fig. 3 we have plotted the differential cross section as a function of the transverse 
momentum of the charged lepton plj, and its pseudorapidity rf , generated with the procedure 
described above, in solid red lines. For comparison we have plotted, in dashed blue lines, the 
NLO differential cross section for W~bb — > l~vhh production of ref. [14], as implemented 
in MCFM 6.0 [48], where spin correlations of the decay products have been included exactly. 
The branching ratio applied to the POWHEG hardest-emission results has been taken from 
the ratio 

BR{W^lv)= = 0.103. (3.22) 

ClO, POWHEG 

As can be seen from the insert in the lower part of the two plots, where we have plotted 
the ratio of the POWHEG hardest emission results over the exact NLO ones (PW/NLO), the 
approximated decay distributions are in very good agreement with the exact ones. 
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Figure 3: Transverse momentum, pip, and pseudorapidity, 77 , distributions for tlie final-state 
lepton produced in W^bb ^ l^D bb, at the LHC with ^/s = f4 TeV. The distributions generated by 
the PDWHEG BOX according to the method described in sec. 3.5 are shown in solid red lines, while 
the NLO results computed with the MCFM code are shown in dashed blue lines. The lower inserts 
show the ratio between the two distributions: POWHEG/NLO. 

4. Phenomenology 

In this section we present a few results for W~bb — t- l~vbh production at the Tevatron and 
at the LHC. Similar results can be obtained for W^bh — t- l^ubb. 

For future reference, we list here the values of all the parameters and physical quantities 
that enter the calculation: 

mw = 80.41 GeV , = 4.62 GeV , mt = 173.1 GeV , 

TvK = 2.141 GeV, BR(Ty ^ /z^) = 0.103 , (4.1) 

and 

sin^ ew = 0.223 , a = 1/132.088832 , Gf = 1.16639 x 10"^ GeV"^ , (4.2) 

from which we derive [g^r = Sm^rGp / V^) 

gyy = 0.6532 . (4.3) 

We have used the CTEQ6.6 pdf set [49], and we have set the renormalization and factor- 
ization scale to the fixed value 

/i = mvi/ + 2 TTife , (4.4) 

from which we compute the two-loop MS strong coupling constant as{fi) = 0.1183 with 5 
light flavors. The VF-boson couplings to quarks are proportional to the Cabibbo-Kobayashi- 
Maskawa (CKM) matrix elements. We use non-zero CKM matrix elements for the first 
two quark generations, Vud = Vcs = 0.974 and Vus = Vcd = 0.227, while we neglect the 
contribution of the third generation, since it is suppressed either by the initial-state quark 
pdfs or by the corresponding CKM matrix elements. 
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Jets are defined using the anti-A;^ algorithm [46] with R = 0.4 and A;Tmin = 5 GeV, 
and are recombined using the default E scheme. 

Since there are no data for W plus two b jets to compare our predictions with, and 
consequently no experimental analysis is available, we have chosen a set of cuts that, while 
reasonable, are less stringent than the experimental ones. In fact, the purpose of this section 
is to show the differences between several POWHEG results, obtained with different showering 
programs, rather then provide predictions for experimentalists, who can use the PDWHEG 
BOX by themselves to generate events and analyze them according to their experimental 
selection criteria. 

The set of cuts for the Tevatron, = 1.96 TeV, and the LHC, a/s = 14 TeV, are the 
following: 

> 15 GeV , \r]^\ < 3 , ^4 > 15 GeV , < 3 , 

p[ > 15 GeV , |r/'| < 3 , I^t > 15 GeV . (4.5) 

We keep only events with at least two 6-jets that pass the cuts on the transverse momentum 
and on the pseudorapidity rj^, disregarding all the other events. Non-6 jets are required 
to have a minimum transverse momentum of 15 GeV and to be in the rapidity region 
\y-^ < 3. Since we have decayed the W boson, we have cuts on the transverse momentum 
pij, and pseudorapidity of the charged lepton, and a cut on the missing energy due 
to the presence of the undetected neutrino. 

In the following, we present several kinematic distributions both for the Tevatron and 
the LHC, and we plot results for the hardest-emission POWHEG cross section with no shower 
(dotted black lines), and for the same results showered by PYTHIA (red solid curves) and 
HERWIG (dashed blue lines). We have run PYTHIA with the Perugia tuning, switching off 
multi-particle interactions'^, in order to make a fair comparison with HERWIG, that uses the 
separate package JIMMY [50] to generate multi-particle interactions. We have run HERWIG 
in its default configuration, with intrinsic px-spreading of 2.5 GeV. 

Our analysis is based on a sample of 29 million events for the Tevatron and 26 million 
events for the LHC, generated with the POWHEG BOX with no folding (see ref. [36] for more 
details). For the Tevatron we got a 15% fraction of negative- weight events and 12% for 
the LHC, that we have kept in our analysis. If one is interested in the generation of only 
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Table 2: Negative- weight fractions as a function of the folding in the radiation variables $rad = 
{^,y,(i)}, at the LHC. 

positive-weight events, then the foldings of the radiation variables ($rad = {C; (p}) should 
^We have set mstp(81)=20 in the code, after the call to pytune. 
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be increased. We have collected in table 2 a few results for different values of the foldings 
in the ^, y and variables. From the table it is evident that, as the product of the folding 
numbers increases, the fraction of negative-weight events gets smaller and smaller, but at 
the price of a decrease in the speed of the code. 



4.1 Tevatron results 




Figure 4: Differential cross sections as a function of the transverse momentum of the liardest b jet, 
leading, the pseudorapidity of the second hardest b jet, 77'' subleading, the invariant mass of the 
leading and subleading b jets, m(,h, and their angular distance Rbb, for W^bb — > l~Dbb production 
at the Tevatron. The different curves represent the results for the PDWHEG hardest emission (dotted 
black), and for POWHEG interfaced with either PYTHIA (sohd red) or HERWIG (dashed blue). 

In fig. 4 we have plotted the differential cross sections as a function of the transverse 
momentum of the b jet with the hardest Pt (called leading b jet), of the pseudorapidity of 
the b jet with the second hardest transverse momentum (called subleading b jet), of the 
invariant mass of the leading and subleading b jets and of the angular distribution R^fj, 
defined as 

Rbb = ^Ayl + AcPl,, (4.6) 

where Ayhh and A0f,b are the difference in the rapidities and in the azimuthal angles of the 
two highest-px b jets. Error bars from the Monte Carlo generator are shown too. 
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Figure 5: Differential cross sections as a function of the transverse momentum, p^, and the 
rapidity, , of the hardest radiated non-6 jet, for W~bb — > l~Dbb production at the Tevatron. The 
different curves represent the resuhs for the POWHEG hardest emission (dotted black), and for POWHEG 
interfaced with either PYTHIA (solid red) or HERWIG (dashed blue). 




Figure 6: Differential cross sections as a function of the transverse momentum, p^, and the 
pseudorapidity, 77', of the lepton for W~bb — ^ l^vbb production at the Tevatron. The different curves 
represent the results for the POWHEG hardest emission (dotted black), and for POWHEG interfaced with 
either PYTHIA (solid red) or HERWIG (dashed blue). 



In fig. 5, we have plotted the differential cross sections as a function of the transverse 
momentum of the hardest radiated (non-6) jet, p^) ai^d its rapidity , while in fig. 6 we have 
plotted the cross sections as a function of the transverse momentum pl^ and pseudorapidity 
rf of the hardest charged lepton. Finally, in fig. 7, we show the differential cross sections 
as a function of {jf^^^ — y-'), i-e. the difference between the rapidity of the Wbb system (the 
W momentum is taken from the showering program, since the neutrino goes undetected) 
and the rapidity of the hardest jet, and as a function of the difference in the azimuthal 
angles of the two b jets, Ac/);,;,. 

In all the plots, we can see that the different showers implemented by PYTHIA and 
HERWIG give slightly different distributions, and we can consider this difference as a the- 
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Figure 7: Differential cross sections as a function of the rapidity difference between the Wbb 
system and the hardest radiated non-6 jet, {y^^^ ~ U''): and the azimuthal angle difference between 
the two b jets A0f,f,, for W^bb — >■ l~Dbb production at the Tevatron. The different curves represent 
the results for the POWHEG hardest emission (dotted black), and for POWHEG interfaced with either 
PYTHIA (sohd red) or HERWIG (dashed blue). 



oretical error associated to showering effects. The trend of the distributions is the same 
in all the plots: the differential cross sections from POWHEG followed by the shower done 
by PYTHIA (POWHEG+PYTHIA) are slightly larger than the corresponding ones showered by 
HERWIG (POWHEG-I- HERWIG). A consequence of this fact is that the cross sections, after the 
cuts of eq. (4.5), are given by 

CpOWHEG = 0.335 pb, 0"powHEG+PYTHIA = 0.291 pb, 0"powHEG+HERWIG = 0.262 pb. (4.7) 

In order to evidentiate the differences between the POWHEG (PW), the POWHEG+PYTHIA 
(PW+PY) and the POWHEG+HERWIG (PW+HW) results, in fig. 8 we have plotted their ratio 
for a sample of distributions. While the ratios of the POWHEG hardest emission cross sec- 
tions over the POWHEG showered results (PW/ (PW+PY) for the shower done by PYTHIA, in 
dashed red lines, and PW/ (PW+HW) for the shower done by HERWIG, in dotted blue lines) 
are just an indication of the effect of the completion of the shower on the POWHEG hardest- 
emission events, the ratios (PW+PY) / (PW+HW) carry information on the effects of the two 
different showering algorithms implemented in PYTHIA and HERWIG. As can be inferred from 
the (PW+PY) / (PW+HW) ratios in the figures, these effects amount to differences of the order 
of 10-20%, and this can be taken as the theoretical errors connected with using different 
showering programs. In general, the ratios are almost flat. The only distribution that 
shows some phase-space dependence is the rapidity of the hardest jet, in the lower left plot 
of fig. 8, corresponding to the ratios of the curves in the right-hand-side plot in fig. 5. Here, 
PW+HW jets tend to be more central in rapidity then the PW+PY ones. We expect this same 
behaviour to be present in fig. 7 as well, since we are plotting the differential distribution 
as a function of [y^^'' — y-?), and the rapidity of the jet enters directly in this quantity. 

4.2 LHC results 

In figs. 9, 10, 11 and 12, we have plotted the same differential cross sections we have 
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Figure 8: Ratios of the differential cross sections for W bb ^ I vbb production at the Tevatron: 
(POWHEG+PYTHIA)/ (POWHEG+HERWIG) in sohd black hnes, POWHEG/(PDWHEG+PYTHIA) in dashed red 
lines and POWHEG/ (POWHEG+HERWIG) in dotted blue hnes. Starting from the upper left corner and 
moving clockwise we show the ratio of the differential cross sections as function of the transverse 
momentum of the hardest b jet, leading, the pseudorapidity of the second hardest b jet, rf' 
subleading, the azimuthal angular difference between the two b jets, Acjibb, and the rapidity of the 
hardest radiated non 6 jet, . 

studied for the Tevatron, this time for the LHC with y/s = 14 TeV. The behaviour of all 
the distributions is the same as for the Tevatron, and the corresponding cross sections after 
cuts are given by 

CpQWHEG = 3.08 pb, CTPOWHEG+PYTHIA = 2.83 pb, TpQWHEG+HERWIG = 2.68 pb. (4.8) 

In fig. 13, we have plotted the ratios (PW+PY) / (PW+HW) in sohd black lines, PW/ (PW+PY) 
in dashed red lines and PW/ (PW+HW) in dotted blue lines. The effects of the two different 
showers, i.e. the ratios (PW+PY) / (PW+HW) are of the order of less than 10% for most of 
the distributions considered, so that the differences between the two showering algorithms 
is less pronounced at the LHC than at the Tevatron. Again, the distribution that turns 
out to be more sensitive to the showering procedure is the rapidity of the hardest jet, as 
illustrated in the lower left-hand-side plot of the figure, where jets from the HERWIG shower 
tend to be more central in rapidity than jets from the PYTHIA shower. 
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Figure 9: Differential cross section as a function of the transverse momentum of the hardest b jet, 
leading, the pseudorapidity of the second hardest b jet, rj^ subleading, the invariant mass of the 
leading and subleading 6 jets rribh, and their angular distance Rbb, for W~bb — l~vbb production at 
the LHC with ^/s = 14 TeV. The different curves represent the results of POWHEG hardest emission 
(dotted black), and of PQWHEG interfaced with either PYTHIA (sohd red) or HERWIG (dashed blue). 

5. Conclusions 

In this article we have presented a next-to-leading order plus parton shower simulation of 
the production of a boson in association with a massive 66 pair, based on the POWHEG 
formalism, with the leptonic decay of the W boson taken into account using standard 
approximated techniques. We have assembled our generator with the aid of the POWHEG 
BOX toolkit [36]. The NLO virtual corrections have been taken from [11, 12, 13] and 
their validity has been expanded in order to account for the case of an off-shell W boson 
production. 

We have validated the code after taking care of activating the POWHEG BOX mechanism 
to protect form Born-zero configurations, and after separating out from the part of the 
real-radiation contribution treated with Monte Carlo techniques the region of hard gluons 
collinear to the final-state massive 6 quarks, to be treated with standard NLO techniques. 
This was done to prevent the enhancement of mass logarithmic terms by the POWHEG factor 
B/B. 
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Figure 10: Differential cross section as a function of the transverse momentum, pip, and the 
rapidity of the hardest radiated non-6 jet, for W^hh l^vbh production at the LHC with 
^/s — 14 TeV. The different curves represent the results of POWHEG hardest emission (dotted black), 
and of POWHEG interfaced with either PYTHIA (sohd red) or HERWIG (dashed blue). 




Figure 11: Differential cross section as a function of the transverse momentum, p^, and the 
pseudorapidity rf of the lepton for W~hh — l~vbb production at the LHC with -y/s — 14 TeV. 
The different curves represent the results of POWHEG hardest emission (dotted black), and of POWHEG 
interfaced with either PYTHIA (solid red) or HERWIG (dashed blue). 



Finally, we have showered the hardest-emission results generated by POWHEG with two 
popular shower Monte Carlo programs: PYTHIA and HERWIG. Looking at various kinematic 
distributions, we have found discrepancies of the order of 10-20% for the Tevatron and of 
less than 10% for the LHC between the two shower Monte Carlo programs. Discrepancies 
larger than the quoted values can be found in distributions involving the rapidity of the 
hardest radiated jet. These discrepancies can be considered as theoretical errors associated 
with the two different showering algorithms. 

The tool we provide will be very important for both Higgs-boson and beyond the 
Standard Model searches at both the Tevatron and the LHC. Indeed W + b jets is one of 
the main backgrounds to these searches and Wbb production is the main contribution. 
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Figure 12: Differential cross section as a function of the rapidity difference between the Wbb 
system and the hardest radiated non-6 jet, {y^^^ — y-'), and the azimuthal angle difference between 
the two b jets A^b^, for W~bb — > l^vbb production at the LHC with -^s — 14 TeV. The different 
curves represent the results of POWHEG hardest emission (dotted black), and of POWHEG interfaced 
with either PYTHIA (sohd red) or HERWIG (dashed blue). 

The code of our generator can be accessed in the POWHEG BOX svn repository 

svn : / / powhegbox . mib . inf n . it/trunk/POWHEG-BOX, 
with username anonymous and password anonymous. 
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